function dataset = scorePeaks(dataset, times, nrOfEpochs, cfg)

%% Baseline subtraction
baselineTf = (times>-4000 & times<-2000) | (times>2000 & times<4000);
M = mean(dataset.mean(baselineTf));
% if M<0.9999 || M>1.0001
%     error ('Baseline check failed'); % Maartje 5 aug 2010
% end;
if exist('cfg', 'var') && cfg.SubtractBaselineBefore30Percent
    m = dataset.mean - M;
else
    m = dataset.mean - 1;
end

%% Smoothing
if exist('cfg', 'var') && strcmp(cfg.Smoothing.Type, 'Gaussian')
    % Smooth the signal before peakfinding.. (DEFAULT is 'none')
    m = GaussianSmooth(m, cfg.Smoothing.Width);
end

%% Scoring
[pks,locsPos] = findpeaks(m, 'minpeakdistance', 100);
[pks,locsNeg] = findpeaks(-m, 'minpeakdistance', 100);
[locs, shuffle] = sort([locsPos, locsNeg]);
type = char(ones(size(locs))*'0');
type(shuffle<=length(locsPos)) = '+';
type(shuffle>length(locsPos)) = '-';
dataset.peaklocs = locs;
dataset.peaktype = type;

%% Calculating significance
if exist('nrOfEpochs', 'var')
    p = [];
    for r = 1:length(dataset.mean)
        z(r) = (abs(dataset.mean(r)-1))/(dataset.std(r)/sqrt(str2double(nrOfEpochs))); % calculate z-value for each latency, based on standard deviation over trials
        p(r,1) = 1 - normcdf(z(r), 0, 1); % calculate p-value
    end;
    dataset.pvalue = p(locs);
end

%% 30%-of-peak-value
%         colr = 'b';
%         plot (times, m, colr);
%         set (gca, 'xtick', floor(times(locs)));
%         grid on;
%         hold on;
dataset.firstEffect = [];
for k = 1:length(locs)
    if type(k) == '-'
        threshTf = m(1:locs(k)) < 0.3*m(locs(k));
    elseif type(k) == '+'
        threshTf = m(1:locs(k)) > 0.3*m(locs(k));
    end
    loc_30perc = find(~threshTf, 1, 'last');
    if isempty(loc_30perc);
        dataset.firstEffect(k).percent30_location = NaN;
        dataset.firstEffect(k).percent30_latency = NaN;
        dataset.firstEffect(k).percent30_peaksSkipped = NaN;
    else
        dataset.firstEffect(k).percent30_location = loc_30perc;
        dataset.firstEffect(k).percent30_latency = times(loc_30perc);
        dataset.firstEffect(k).percent30_peaksSkipped = sum(locs>loc_30perc & locs<locs(k));
    end;
    %             plot(times([loc_30perc loc_30perc]), [min(m) max(m)], colr);
end
end